;; make initial data for  Jablonowski & Williamson (2006) dynamical core test

load "$NCARG_ROOT/lib/ncarg/nclscripts/csm/contributed.ncl"

fout_steady="./ini_steady.t30.grd"
fnc_steady =addfile("./ini_steady.t30.nc","c")
fout_barocl="./ini_barocl.t30.grd"
fnc_barocl =addfile("./ini_barocl.t30.nc","c")


nlat=48
nlon=96
nlev=7

;; set coordinates
gau_info=gaus(nlat/2)
lat=tofloat(gau_info(:,0))
dlon=360./nlon
lon=dlon*ispan(0,nlon-1,1)
sig=(/ 0.950, 0.835, 0.685, 0.510, 0.340, 0.200, 0.080 /)

;; set-up variables
u = new((/nlev,nlat,nlon/),"float")
u!0="lev"
u!1="lat"
u!2="lon"
u&lev=sig
u&lat=lat
u&lon=lon
v = u
t = u
q = u

;; surface pressure
ps = new((/nlat,nlon/),"float")
ps!0="lat"
ps!1="lon"
ps&lat=lat
ps&lon=lon
;; rain (dummy)
rain=ps

;; parameters
Omega=7.29212e-5
a=6.371229e+6
g=9.80616
Rd=287.
Gamma=0.005

u0=35.
T0=288.
dT=4.8e+5

pi=4*atan2(1,1)
eta_0=0.252
eta_t=0.2
eta_v=(sig-eta_0)*(pi/2)  ;; Eq.(1)
c = cos(eta_v)^1.5
phi=lat*(pi/180)
lmd=lon*(2.*pi/360.)
cosphi=cos(phi)
sinphi=sin(phi)


eta_v3d = conform(u, eta_v, 0)
c3d     = conform(u,     c, 0)
sig3d   = u
sig3d   = conform(u,   sig, 0)
cosphi3d= conform(u,cosphi, 1)
sinphi3d= conform(u,sinphi, 1)
lmd3d   = conform(u,   lmd, 2)

;; zonal wind
u = u0*c3d*(4*(sinphi3d^2)*(cosphi3d^2))  ;; Eq.(2)

;; meridional wind
v = 0

;; temperature
t=T0*(sig3d^(Rd*Gamma/g))  ;; Eq.(4)
t({eta_t:0},:,:) = t({eta_t:0},:,:) + dT*((eta_t-sig3d({eta_t:0},:,:))^5) ;; Eq.(5)
;; Eq.(6)
t=t+(3./4.)*(sig3d*pi*u0/Rd)*sin(eta_v3d)*sqrt(cos(eta_v3d))\
           *( (-2*(sinphi3d^6)*((cosphi3d^2)+1./3.)+10./63.)*2*u0*c3d\
               + ((8./5.)*(cosphi3d^3)*((sinphi3d^2)+2./3.)-pi/4.)*a*Omega )

;; humidity
q=0
;; surface pressure
ps = 1e+5
;; rain
rain=0

;;; output  steady state to file
setfileoption("bin","WriteByteOrder","BigEndian")
fbindirwrite(fout_steady,   u)
fbindirwrite(fout_steady,   v)
fbindirwrite(fout_steady,   t)
fbindirwrite(fout_steady,   q)
fbindirwrite(fout_steady,  ps)
fbindirwrite(fout_steady,rain)
fnc_steady->u=u
fnc_steady->t=t

;;;;;;;;; baroclinic wave ;;;;;;;;
lmd_c=pi/9.
phi_c=(2.*pi)/9.
up=1.
R=a/10.

r = a*acos( sin(phi_c)*sinphi3d + cos(phi_c)*cosphi3d*cos(lmd3d-lmd_c) ) ;; Eq.(11)
u = u + up*exp(-((r/R)^2)) ;; Eq.(10)

;;; output  baroclinic initial condition to file
setfileoption("bin","WriteByteOrder","BigEndian")
fbindirwrite(fout_barocl,   u)
fbindirwrite(fout_barocl,   v)
fbindirwrite(fout_barocl,   t)
fbindirwrite(fout_barocl,   q)
fbindirwrite(fout_barocl,  ps)
fbindirwrite(fout_barocl,rain)
fnc_barocl->u=u
fnc_barocl->t=t
